{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.interpolate import Rbf\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "np.random.seed(1981)\n",
    "\n",
    "x, y, z = np.random.random((3, 10))\n",
    "xi, yi = np.mgrid[0:1:100j, 0:1:100j]\n",
    "func = Rbf(x, y, z, function='linear')\n",
    "zi = -1*func(xi, yi)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# Plot flowlines\n",
    "dy, dx = np.gradient(-zi.T)\n",
    "\n",
    "# Contour gridded head observations\n",
    "contours = ax.contour(xi, yi, zi[:,::-1], linewidths=2)\n",
    "ax.clabel(contours)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_idx(np_linarray, val):\n",
    "    i = 0\n",
    "    if (val < np_linarray[0]) or (val > np_linarray[-1]):\n",
    "        i = -1\n",
    "    else:\n",
    "        for i in range(len(np_linarray)-1):\n",
    "            if val < np_linarray[i+1] and val > np_linarray[i]:\n",
    "                break\n",
    "        \n",
    "    return i\n",
    "\n",
    "def get_grad(X, Y, z, pos):\n",
    "    dy, dx = np.gradient(-z.T)\n",
    "    grad = [np.nan, np.nan]\n",
    "    \n",
    "    idx_x = get_idx(X[:,0], pos[0])\n",
    "    idx_y = get_idx(Y[0,:], pos[1])    \n",
    "    grad = np.array([ dx[idx_y,idx_x] , dy[idx_y,idx_x] ])\n",
    "    #print(idx_x, idx_y, pos[0], pos[1], grad, zi[idx_x, idx_y])\n",
    "    return grad\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "position_init = np.array([0.7, 0.6]) # other init position?\n",
    "velocity = np.array([-0.01, -0.01])\n",
    "alpha = 0.5\n",
    "learning_rate = 0.05                 # high or low learning rate?\n",
    "eproch = 10000                       # less eproch?\n",
    "delta = 1e-7\n",
    "rho1 = 0.9\n",
    "rho2 = 0.999\n",
    "\n",
    "def func_sgd(theta, eproch, xi, yi, zi, learning_rate):\n",
    "    l_posx = []\n",
    "    l_posy = []\n",
    "    for i in range(eproch):\n",
    "        grad = get_grad(xi, yi, zi, theta)\n",
    "        theta = theta-grad*learning_rate\n",
    "        l_posx.append(theta[0])\n",
    "        l_posy.append(1-theta[1])\n",
    "    \n",
    "    return l_posx, l_posy\n",
    "\n",
    "def func_momentum(theta, eproch, xi, yi, zi, learning_rate, alpha, velocity):\n",
    "    l_posx = []\n",
    "    l_posy = []\n",
    "    for i in range(eproch):\n",
    "        grad = get_grad(xi, yi, zi, theta)\n",
    "        velocity = alpha*velocity - grad*learning_rate\n",
    "        theta  = theta + velocity\n",
    "        l_posx.append(theta[0])\n",
    "        l_posy.append(1-theta[1])\n",
    "\n",
    "    return l_posx, l_posy\n",
    "\n",
    "\n",
    "def func_adagrad(theta, eproch, xi, yi, zi, learning_rate, delta):\n",
    "    l_posx = []\n",
    "    l_posy = []\n",
    "    r = 0\n",
    "    for i in range(eproch):\n",
    "        grad = get_grad(xi, yi, zi, theta)\n",
    "        r += np.dot(grad, grad)\n",
    "        theta  = theta- learning_rate/(delta+np.sqrt(r)) * grad\n",
    "        l_posx.append(theta[0])\n",
    "        l_posy.append(1-theta[1])\n",
    "\n",
    "    return l_posx, l_posy\n",
    "\n",
    "\n",
    "def func_adam(theta, eproch, xi, yi, zi, learning_rate, delta, rho1, rho2):\n",
    "    l_posx = []\n",
    "    l_posy = []\n",
    "    s = 0\n",
    "    r = 0\n",
    "    for i in range(eproch):\n",
    "        grad = get_grad(xi, yi, zi, theta)\n",
    "        t = i+1\n",
    "        s = rho1 * s + (1-rho1) * grad\n",
    "        r = rho2 * 2 + (1-rho2) * np.dot(grad, grad)\n",
    "        s_hat = s / (1-rho1**t)\n",
    "        r_hat = r / (1-rho2**t)\n",
    "        \n",
    "        theta  = theta- learning_rate/(delta+np.sqrt(r_hat)) * s_hat\n",
    "        l_posx.append(theta[0])\n",
    "        l_posy.append(1-theta[1])\n",
    "\n",
    "    return l_posx, l_posy\n",
    "\n",
    "\n",
    "l_posx_gd,l_posy_gd = func_sgd(\n",
    "                        position_init, eproch, \n",
    "                        xi, yi, zi, learning_rate)\n",
    "\n",
    "l_posx_momentum,l_posy_momentum = func_momentum(\n",
    "                        position_init, eproch, \n",
    "                        xi, yi, zi, learning_rate, alpha, velocity)\n",
    "\n",
    "l_posx_adagrad,l_posy_adagrad = func_adagrad(\n",
    "                        position_init, eproch, \n",
    "                        xi, yi, zi, learning_rate, delta)\n",
    "\n",
    "l_posx_adam,l_posy_adam = func_adam(\n",
    "                        position_init, eproch, \n",
    "                        xi, yi, zi, learning_rate, delta, rho1, rho2)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "contours = ax.contour(xi, yi, zi[:,::-1], linewidths=2)\n",
    "ax.clabel(contours)\n",
    "ax.plot(l_posx_gd, l_posy_gd, \".\", label=\"SGD\")\n",
    "ax.plot(l_posx_momentum, l_posy_momentum, \"g.\", label=\"Momentum\")\n",
    "ax.plot(l_posx_adagrad, l_posy_adagrad, \"y.\", label=\"Adagrad\")\n",
    "ax.plot(l_posx_adam,    l_posy_adam,    \"c.\", label=\"Adam\")\n",
    "\n",
    "ax.plot(l_posx_gd[-1],       l_posy_gd[-1],       \"ro\")\n",
    "ax.plot(l_posx_momentum[-1], l_posy_momentum[-1], \"ro\")\n",
    "ax.plot(l_posx_adagrad[-1],  l_posy_adagrad[-1],  \"ro\")\n",
    "ax.plot(l_posx_adam[-1],     l_posy_adam[-1],     \"ro\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
